1 CONFIGURATIONS

# Set the working directory and tool paths on your local computer.
WD <- '/Users/Alec/motrpac/20210826_pass1c-umich'
# Set the gsutil path

knitr::opts_knit$set(root.dir=WD)
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(warning = FALSE)
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(cache = FALSE)

1.0.0.0.1 CSS Top Styling

write_css = FALSE
if(write_css){
  writeLines("td, th { padding : 3px } th { background-color:white; color:black; border:1px solid black; text-align:center } td {color:black; border:1px solid black; word-wrap:break-word; white-space:nowrap; overflow: hidden; text-overflow: ellipsis; max-width:300px; text-align:left}", con=file.path(normalizePath(WD), "style/style.css"))
}

1.0.0.0.2 Overview

1.0.0.0.3 Goals of Analysis

Examine data dimensions
Examine zero, negative, and missing Values
Impute missing values
Visualize sample-by-sample correlations
Identify outliers
Identify major causes of variance (drift, sex, control group, hour)
Visualize sample-by-feature heatmaps under different transformations
Examine sample median and sd distributions between each transformation
Transform data
Remove outlier samples
Save the processed data

2 Prepare Environment & Load Data

2.1 Setup the Environment

################################################################################
##### Resources and Dependencies ###############################################
################################################################################
# Whether to knit document and display data
knit_time = TRUE

# Load dependencies
pacs...man <- c("tidyverse","kableExtra","devtools","MotrpacBicQC","impute","glue",
                "rethinking")
for(pac in pacs...man){
  suppressWarnings(suppressPackageStartupMessages(library(pac, character.only = TRUE)))
}

#browseVignettes("MotrpacBicQC")
############################################################
##### Functions ############################################
############################################################

# Name functions
select <- dplyr::select
map <- purrr::map
desc <- dplyr::desc
arrange <- dplyr::arrange
melt <- reshape2::melt
mutate <- dplyr::mutate
glue <- glue::glue

# Global options
options(dplyr.print_max = 100)
options(scipen=10000)

# Colors
redblue100<-rgb(read.table(paste0(WD,'/colors/redblue100.txt'),sep='\t',row.names=1,header=T))

2.2 Log Variables (1)

ds <- 'pass1a'
site <- 'umich'
tech <- 'rpneg'
TIS <- list('PLA', 'HIP', 'GAS', 'HRT', 'KID', 'LUN', 'LIV', 'BADI', 'WADI')
tis <- list('pla', 'hip', 'gas', 'hrt', 'kid', 'lun', 'liv', 'badi', 'wadi')

2.3 Load Phenotype data

# Phenotype Data (1A)
#########################
pheno_file <- glue("{WD}/data/20201021_pass1a-06-pheno-viallabel_steep.txt")
pheno_df1a <- read.csv(pheno_file, header = T, sep = '\t')
pheno_df1a <- pheno_df1a %>%
  mutate(tissue = case_when(Specimen.Processing.sampletypedescription == 'Brown Adipose' ~ 'badi',
                            Specimen.Processing.sampletypedescription == 'EDTA Plasma' ~ 'pla',
                            Specimen.Processing.sampletypedescription == 'Gastrocnemius' ~ 'gas',
                            Specimen.Processing.sampletypedescription == 'Heart' ~ 'hrt',
                            Specimen.Processing.sampletypedescription == 'Kidney' ~ 'kid',
                            Specimen.Processing.sampletypedescription == 'Liver' ~ 'liv',
                            Specimen.Processing.sampletypedescription == 'Lung' ~ 'lun',
                            Specimen.Processing.sampletypedescription == 'White Adipose' ~ 'wadi',
                            Specimen.Processing.sampletypedescription == 'PaxGene RNA' ~ 'pax',
                            Specimen.Processing.sampletypedescription == 'Hippocampus' ~ 'hip',
                            Specimen.Processing.sampletypedescription == 'Cortex' ~ 'cor',
                            Specimen.Processing.sampletypedescription == 'Hypothalamus' ~ 'hyp',
                            Specimen.Processing.sampletypedescription == 'Vastus Lateralis' ~ 'vas',
                            Specimen.Processing.sampletypedescription == 'Tibia' ~ 'tib',
                            Specimen.Processing.sampletypedescription == 'Aorta' ~ 'aor',
                            Specimen.Processing.sampletypedescription == 'Small Intestine' ~ 'sma',
                            Specimen.Processing.sampletypedescription == 'Adrenals' ~ 'adr',
                            Specimen.Processing.sampletypedescription == 'Colon' ~ 'col',
                            Specimen.Processing.sampletypedescription == 'Spleen' ~ 'spl',
                            Specimen.Processing.sampletypedescription == 'Testes' ~ 'tes',
                            Specimen.Processing.sampletypedescription == 'Ovaries' ~ 'ova'))
pheno_df1a$viallabel <- as.character(pheno_df1a$viallabel)

2.4 Load Data

2.4.1 Load the sample order files

# #created in 20210910_pass1a-umich-sample-annotation_steep.Rmd
# order_file <- glue("{WD}/data/20200910_pass1a1c-sample-order_steep.txt")
# sample.order<-read.delim(order_file,header=T, sep="\t")

# Load the prior pass1a data (takes a few minutes)
pass1a_nested_file <- glue("{WD}/../20200915_metabolomics-pass1a/data/20201010_pass1a-metabolomics-countdata-nested_steep.rds")
pass1a_df <- readRDS(pass1a_nested_file)

# pla
pla_rpneg_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rpneg') %>%
  filter(TISSUE == 'plasma') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)

# hip
hip_rpneg_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rpneg') %>%
  filter(TISSUE == 'hippocampus') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)

# gas
gas_rpneg_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rpneg') %>%
  filter(TISSUE == 'gastrocnemius') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)

# hrt
hrt_rpneg_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rpneg') %>%
  filter(TISSUE == 'heart') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)

# kid
kid_rpneg_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rpneg') %>%
  filter(TISSUE == 'kidney') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)

# lun
lun_rpneg_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rpneg') %>%
  filter(TISSUE == 'lung') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)

# liv
liv_rpneg_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rpneg') %>%
  filter(TISSUE == 'liver') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)

# badi
badi_rpneg_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rpneg') %>%
  filter(TISSUE == 'brown-adipose') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)

# wadi
wadi_rpneg_pass1a.0.order <- pass1a_df %>%
  filter(STUDY_INSTITUTE == 'University of Michigan') %>%
  filter(NAMED == 'named') %>%
  filter(DATASET == 'untargeted') %>%
  filter(METAB_FAMILY == 'rpneg') %>%
  filter(TISSUE == 'white-adipose') %>%
  select(SAMPLE_DATA) %>%
  unnest(SAMPLE_DATA) %>%
  left_join(y = pheno_df1a, by = c('sample_id' = 'viallabel')) %>%
  select(sample_id, sample_type, sample_order, Key.anirandgroup, Registration.sex)

rm(pass1a_df)

# Create a list of annotation matrixes
rpneg_pass1a.0.order <- list(pla_rpneg_pass1a.0.order, hip_rpneg_pass1a.0.order, gas_rpneg_pass1a.0.order, 
                 hrt_rpneg_pass1a.0.order, kid_rpneg_pass1a.0.order, lun_rpneg_pass1a.0.order, 
                 liv_rpneg_pass1a.0.order, badi_rpneg_pass1a.0.order, wadi_rpneg_pass1a.0.order)

# Save the Sample order file
save(pla_rpneg_pass1a.0.order, hip_rpneg_pass1a.0.order, gas_rpneg_pass1a.0.order, 
                 hrt_rpneg_pass1a.0.order, kid_rpneg_pass1a.0.order, lun_rpneg_pass1a.0.order, 
                 liv_rpneg_pass1a.0.order, badi_rpneg_pass1a.0.order, wadi_rpneg_pass1a.0.order,
     file = glue("{WD}/data/UM-rpneg/rpneg_pass1a.0.order.RData"))

2.4.2 Load the matrices as .RData files

# UM-rpneg
load(file = glue("{WD}/data/UM-rpneg/UM_rpneg.0.RData"))
pla1a.0 <- pla_rpneg_pass1a.0
hip1a.0 <- hip_rpneg_pass1a.0
gas1a.0 <- gas_rpneg_pass1a.0
hrt1a.0 <- hrt_rpneg_pass1a.0
kid1a.0 <- kid_rpneg_pass1a.0
lun1a.0 <- lun_rpneg_pass1a.0
liv1a.0 <- liv_rpneg_pass1a.0
badi1a.0 <- badi_rpneg_pass1a.0
wadi1a.0 <- wadi_rpneg_pass1a.0

# Create a list of matrices
pass1a.0 <- list(pla1a.0, hip1a.0, gas1a.0, hrt1a.0, kid1a.0, lun1a.0, liv1a.0, badi1a.0, wadi1a.0)

2.5 Remove Internal Standards

is <- list()
for(i in 1:length(pass1a.0)){
  # Remove internal standards
  is[[i]] <- colnames(pass1a.0[[i]])[grepl("istd", colnames(pass1a.0[[i]]), ignore.case = TRUE)]
  # Subset matrix
  pass1a.0[[i]] <- pass1a.0[[i]][, !colnames(pass1a.0[[i]]) %in% is]
}

2.5.1 Create Annotations

i <- 1
tmp.iter1a <- tmp.ref1a <- tmp.sample1a <- pass1a.0.nr <- list()
for(i in 1:length(rpneg_pass1a.0.order)){
  print(tis[[i]])
  # Collect the distances from reference samples
  tmp.iter1a[[i]] <- rpneg_pass1a.0.order[[i]] %>%
    mutate(reference = ifelse(str_sub(sample_id,1,1) == '9', 0, 1)) %>%
    mutate(drift = ifelse(sample_type == 'QC-DriftCorrection', 1, 0)) %>%
    arrange(sample_order)
  tmp.iter1a[[i]]$right_p <- 0
  for(j in 1:nrow(tmp.iter1a[[i]])){
    # set t
    if(tmp.iter1a[[i]][j,'reference'] == 1){
      t = 0
    }else if(tmp.iter1a[[i]][j,'reference'] == 0){
      t = t + 1
    }
    tmp.iter1a[[i]][j,"right_p"] <- t
  }
  t=0
  tmp.iter1a[[i]]$right_p_d <- 0
  for(j in 1:nrow(tmp.iter1a[[i]])){
    # set t
    if(tmp.iter1a[[i]][j,'drift'] == 1){
      t = 0
    }else if(tmp.iter1a[[i]][j,'drift'] == 0){
      t = t + 1
    }
    tmp.iter1a[[i]][j,"right_p_d"] <- t
  }
  t=0

  tmp.iter1a[[i]] <- tmp.iter1a[[i]] %>%
    arrange(desc(sample_order))
  tmp.iter1a[[i]]$left_p <- 0
  for(j in 1:nrow(tmp.iter1a[[i]])){
    # set t
    if(tmp.iter1a[[i]][j,'reference'] == 1){
      t = 0
    }else if(tmp.iter1a[[i]][j,'reference'] == 0){
      t = t + 1
    }
    tmp.iter1a[[i]][j,"left_p"] <- t
  }
  t=0
  tmp.iter1a[[i]]$left_p_d <- 0
  for(j in 1:nrow(tmp.iter1a[[i]])){
    # set t
    if(tmp.iter1a[[i]][j,'drift'] == 1){
      t = 0
    }else if(tmp.iter1a[[i]][j,'drift'] == 0){
      t = t + 1
    }
    tmp.iter1a[[i]][j,"left_p_d"] <- t
  }
  t=0

  tmp.iter1a[[i]] <- tmp.iter1a[[i]] %>%
    rowwise() %>%
    mutate(min_p = min(c(left_p,right_p) ,na.rm = TRUE)) %>%
    mutate(sum_p = sum(c(left_p,right_p) ,na.rm = TRUE)) %>%
    mutate(min_p_d = min(c(left_p_d,right_p_d) ,na.rm = TRUE)) %>%
    mutate(sum_p_d = sum(c(left_p_d,right_p_d) ,na.rm = TRUE)) %>%
    arrange(sample_order)
  tmp.join1a <- tmp.iter1a[[i]] %>%
    select(sample_id, sample_order, left_p, right_p, min_p, sum_p, left_p_d, right_p_d, min_p_d, sum_p_d)

  # Collect the sample order for test+ref samples
  tmp.ref1a[[i]] <- rpneg_pass1a.0.order[[i]] %>%
    arrange(sample_order) %>% 
    mutate(control = ifelse(grepl('Control', Key.anirandgroup), 1, 0)) %>%
    mutate(drift = ifelse(grepl('Drift', sample_type), 1, 0)) %>%
    mutate(reference = ifelse(str_sub(sample_id,1,1) == '9', 0, 1)) %>%
    mutate(time = case_when(grepl('IPE', Key.anirandgroup) ~ 0,
                            grepl('0 hr', Key.anirandgroup) ~ 0,
                            grepl('0.5 hr', Key.anirandgroup) ~ 0.5,
                            grepl('1 hr', Key.anirandgroup) ~ 1,
                            grepl('4 hr', Key.anirandgroup) ~ 4,
                            grepl('7 hr', Key.anirandgroup) ~ 7,
                            grepl('24 hr', Key.anirandgroup) ~ 24,
                            grepl('48 hr', Key.anirandgroup) ~ 48)) %>%
    left_join(y = tmp.join1a) %>%
    select(sample_id,sample_order, Registration.sex, control, time, drift, reference,
         left_p, right_p, min_p, sum_p, left_p_d, right_p_d, min_p_d, sum_p_d)
  N2 <- tmp.ref1a[[i]][,1] %>% unlist()
  miss1 <- N2[!N2 %in% row.names(pass1a.0[[i]])] # Verify all samples are in the pass1a.0[[i]] file
  miss1
  # sample.order %>% filter(sample_id == "90750016606")
  print('Vice Versa:')
  miss2 <- row.names(pass1a.0[[i]])[!row.names(pass1a.0[[i]]) %in% N2]
  miss2
  N2 <- N2[!N2 %in% c(miss1,miss2)] # TODO: investigate why samples are missing, continue for now
  # Reorder pass1a.0[[i]] by run order
  pass1a.0[[i]] <- pass1a.0[[i]][N2,]
  all(N2 == row.names(pass1a.0[[i]])) # Must be true
  tmp.ref1a[[i]] <- tmp.ref1a[[i]] %>%
    filter(sample_id %in% N2) %>%
    select(sample_order, Registration.sex, control, time, drift, reference,
           left_p, right_p, min_p, sum_p, left_p_d, right_p_d, min_p_d, sum_p_d) %>% as.matrix()

  # Collect the sample order for test samples
  options(digits = 14)
  tmp.sample1a[[i]] <- rpneg_pass1a.0.order[[i]] %>%
    filter(substr(sample_id, 1, 1) == '9') %>%
    arrange(sample_order) %>% 
    mutate(control = ifelse(grepl('Control', Key.anirandgroup), 1, 0)) %>%
    mutate(time = case_when(grepl('IPE', Key.anirandgroup) ~ 0,
                            grepl('0 hr', Key.anirandgroup) ~ 0,
                            grepl('0.5 hr', Key.anirandgroup) ~ 0.5,
                            grepl('1 hr', Key.anirandgroup) ~ 1,
                            grepl('4 hr', Key.anirandgroup) ~ 4,
                            grepl('7 hr', Key.anirandgroup) ~ 7,
                            grepl('24 hr', Key.anirandgroup) ~ 24,
                            grepl('48 hr', Key.anirandgroup) ~ 48)) %>%
    left_join(y = tmp.join1a) %>%
    mutate(sample_id = str_replace_all(sample_id, pattern = '-', replacement = '')) %>%
    mutate(sample_id = as.numeric(sample_id)) %>%
    select(sample_id, sample_order, Registration.sex, control, time, 
           left_p, right_p, min_p, sum_p, left_p_d, right_p_d, min_p_d, sum_p_d) %>%
    as.matrix()

  N1 <- row.names(pass1a.0[[i]])[substr(row.names(pass1a.0[[i]]), 1, 1) == '9']
  tmp.sample1a[[i]] <- tmp.sample1a[[i]][tmp.sample1a[[i]][,1] %in% N1,]
  # Reorder pass1a.0[[i]].nr by run order
  pass1a.0.nr[[i]] <- pass1a.0[[i]][N1,]
  tmp.sample1a[[i]][,1][!tmp.sample1a[[i]][,1] %in% row.names(pass1a.0.nr[[i]])] # Verify all samples are in the pass1a.0[[i]] file
  all(as.character(tmp.sample1a[[i]][,1]) == row.names(pass1a.0.nr[[i]]))
  # If out of order, command below will ensure pass1a.0[[i]].nr in run order
  #pass1a.0[[i]].nr <- pass1a.0[[i]].nr[as.character(tmp.sample1a[[i]][,1]),]
}

## [1] "pla"
## [1] "Vice Versa:"
## [1] "hip"
## [1] "Vice Versa:"
## [1] "gas"
## [1] "Vice Versa:"
## [1] "hrt"
## [1] "Vice Versa:"
## [1] "kid"
## [1] "Vice Versa:"
## [1] "lun"
## [1] "Vice Versa:"
## [1] "liv"
## [1] "Vice Versa:"
## [1] "badi"
## [1] "Vice Versa:"
## [1] "wadi"
## [1] "Vice Versa:"

3 Dimensions, Zero/Neg/Missing Values, & Log2

3.1 Dimensions (with reference samples)

# Lists
NR <- P <- list()
for(i in 1:length(pass1a.0)){
  print(tis[[i]])
  NR[[i]] <- dim(pass1a.0[[i]])[1]
  P[[i]] <- dim(pass1a.0[[i]])[2]
  print(dim(pass1a.0[[i]]))
}

## [1] "pla"
## [1] 96 97
## [1] "hip"
## [1] 113 139
## [1] "gas"
## [1] 98 78
## [1] "hrt"
## [1] 120  94
## [1] "kid"
## [1] 120  87
## [1] "lun"
## [1] 120 102
## [1] "liv"
## [1] 103  81
## [1] "badi"
## [1] 116 143
## [1] "wadi"
## [1] 112  78

3.2 Dimensions (without reference samples)

N <- list()
for(i in 1:length(pass1a.0.nr)){
  print(tis[[i]])
  N[[i]] <- dim(pass1a.0.nr[[i]])[1]
  print(dim(pass1a.0.nr[[i]]))
}

## [1] "pla"
## [1] 77 97
## [1] "hip"
## [1]  78 139
## [1] "gas"
## [1] 78 78
## [1] "hrt"
## [1] 78 94
## [1] "kid"
## [1] 78 87
## [1] "lun"
## [1]  78 102
## [1] "liv"
## [1] 78 81
## [1] "badi"
## [1]  78 143
## [1] "wadi"
## [1] 76 78

3.3 Negative or Zero Values

confirmed: no negative or zero values

for(i in 1:length(pass1a.0.nr)){
  print(tis[[i]])
  print(min(pass1a.0[[i]],na.rm=T))
}

## [1] "pla"
## [1] 36.40427817
## [1] "hip"
## [1] 71.37304298
## [1] "gas"
## [1] 49.21625509
## [1] "hrt"
## [1] 120.973489
## [1] "kid"
## [1] 101.3145201
## [1] "lun"
## [1] 61.68128751
## [1] "liv"
## [1] 502.8654934
## [1] "badi"
## [1] 100.201935
## [1] "wadi"
## [1] 63.25290985

3.4 Missing Features (with references)

Blank reference samples at the beginning and end

pass1a.0.f.c0 <- list()
for(i in 1:length(pass1a.0)){
  pass1a.0.f.c0[[i]] <-apply(pass1a.0[[i]],1,function(x) sum(is.na(x))) 
}

# Plot
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  plot(pass1a.0.f.c0[[i]], ylim = c(0,max(unlist(P))), main = glue("{tis[[i]]}"), xlab = 'Samples', ylab = 'Missing Values')
}

3.5 Blank Samples (without references)

No blank test samples

pass1a.0.nr.f.c0 <- list()
for(i in 1:length(pass1a.0.nr)){
  pass1a.0.nr.f.c0[[i]] <-apply(pass1a.0.nr[[i]],1,function(x) sum(is.na(x))) 
}

# Plot
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0.nr)){
  plot(pass1a.0.nr.f.c0[[i]], ylim = c(0,max(unlist(P))), main = glue("{tis[[i]]}"), xlab = 'Samples', ylab = 'Missing Values')
}

3.6 Examine Distribution of missing values (with reference samples)

pass1a.0.f.c0 <- list()
for(i in 1:length(pass1a.0.nr)){
  pass1a.0.f.c0[[i]] <- apply(pass1a.0[[i]],2,function(x) sum(is.na(x))) 
}
# Plot
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0.nr)){
  plot(pass1a.0.f.c0[[i]], ylim = c(0,max(unlist(NR))), main = glue("{tis[[i]]}"), xlab = 'Features', ylab = 'Missing Values')
}

3.7 Examine Distribution of missing values (test samples only)

pass1a.0.nr.f.c0 <- list()
for(i in 1:length(pass1a.0.nr)){
  pass1a.0.nr.f.c0[[i]] <- apply(pass1a.0.nr[[i]],2,function(x) sum(is.na(x))) 
}
# Plot
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0.nr)){
  plot(pass1a.0.nr.f.c0[[i]], ylim = c(0,max(unlist(N))), main = glue("{tis[[i]]}"), xlab = 'Features', ylab = 'Missing Values'); abline(h = 20)
}

3.8 Remove high-missing features

Remove high missing features

rm_n <- list()
for(i in 1:length(pass1a.0.nr.f.c0)){
  print(tis[[i]])
  rm_n[[i]] <- sum(pass1a.0.nr.f.c0[[i]]>=20)
  pass1a.0[[i]] <- pass1a.0[[i]][,pass1a.0.nr.f.c0[[i]]<20]
  pass1a.0.nr[[i]] <- pass1a.0.nr[[i]][,pass1a.0.nr.f.c0[[i]]<20]
  print(dim(pass1a.0[[i]]))
  print(dim(pass1a.0.nr[[i]]))
}

## [1] "pla"
## [1] 96 96
## [1] 77 96
## [1] "hip"
## [1] 113 109
## [1]  78 109
## [1] "gas"
## [1] 98 76
## [1] 78 76
## [1] "hrt"
## [1] 120  94
## [1] 78 94
## [1] "kid"
## [1] 120  87
## [1] 78 87
## [1] "lun"
## [1] 120 102
## [1]  78 102
## [1] "liv"
## [1] 103  81
## [1] 78 81
## [1] "badi"
## [1] 116 123
## [1]  78 123
## [1] "wadi"
## [1] 112  76
## [1] 76 76

3.9 Take the log2

pass1a.0.1 <- pass1a.0.nr1 <- list()
for(i in 1:length(pass1a.0)){
  pass1a.0.1[[i]] <- log(pass1a.0[[i]], 2)
  pass1a.0.nr1[[i]] <- log(pass1a.0.nr[[i]], 2)
}

3.10 Total missing values (includes reference samples)

for(i in 1:length(pass1a.0.1)){
  print(glue("{tis[[i]]}"))
  print(sum(is.na(pass1a.0.1[[i]])))
}

## pla
## [1] 2
## hip
## [1] 408
## gas
## [1] 20
## hrt
## [1] 268
## kid
## [1] 289
## lun
## [1] 363
## liv
## [1] 105
## badi
## [1] 563
## wadi
## [1] 450

3.11 Total missing values (test samples)

feature_impute <- list()
for(i in 1:length(pass1a.0.nr1)){
  print(tis[[i]])
  print(sum(is.na(pass1a.0.nr1[[i]])))
  feature_impute[[i]] <- apply(is.na(pass1a.0.nr1[[i]]),2,sum)[apply(is.na(pass1a.0.nr1[[i]]),2,sum) > 0]
}

## [1] "pla"
## [1] 0
## [1] "hip"
## [1] 3
## [1] "gas"
## [1] 20
## [1] "hrt"
## [1] 19
## [1] "kid"
## [1] 9
## [1] "lun"
## [1] 5
## [1] "liv"
## [1] 80
## [1] "badi"
## [1] 10
## [1] "wadi"
## [1] 45

3.12 Confirm New Distribution of missing features (test samples only)

pass1a.0.nr1.f.c0 <- list()
for(i in 1:length(pass1a.0.nr1)){
  pass1a.0.nr1.f.c0[[i]] <- apply(pass1a.0.nr1[[i]],2,function(x) sum(is.na(x))) 
}

# Plot
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0.nr1)){
  plot(pass1a.0.nr1.f.c0[[i]], ylim = c(0,max(unlist(NR))), main = glue("{tis[[i]]}"), xlab = 'Features', ylab = 'Missing Values'); abline(h = 20)
}

3.13 Log Variables (2)

# NR
# N
# P
neg_vals <- 0
zero_vals <- 0
feature_na_filter <- 20
knn_k <- 10


P1 <- NR1 <- N1 <- na_vals_impute <- list()
for(i in 1:length(pass1a.0)){
  print(tis[[i]])
  P1[[i]] <- dim(pass1a.0.nr[[i]])[2]
  NR1[[i]] <- dim(pass1a.0[[i]])[1]
  N1[[i]] <- dim(pass1a.0.nr[[i]])[1]
  na_vals_impute[[i]] <- sum(is.na(pass1a.0.nr1[[i]]))
  print(feature_impute[[i]])
}

## [1] "pla"
## named integer(0)
## [1] "hip"
## glutamate-13C5 [iSTD] 
##                     3 
## [1] "gas"
## Tetradecadienoic acid      Myristoleic acid         Deoxycytidine 
##                     3                     2                     1 
##     Pentadecylic acid            Xanthosine 
##                     1                    13 
## [1] "hrt"
##        alanine-13C3 [iSTD] 3-Methyl-2-oxovaleric acid 
##                          2                          3 
##                Ketoleucine                Lauric acid 
##                          2                          8 
##            Lignoceric acid       Ursodeoxycholic acid 
##                          3                          1 
## [1] "kid"
##          Lactic acid          Capric acid Pregnenolone sulfate 
##                    2                    4                    2 
##          Cholic acid 
##                    1 
## [1] "lun"
## Tauro-alpha-muricholic acid 
##                           5 
## [1] "liv"
##                      Tetradecadienoic acid 
##                                          2 
##                           Myristoleic acid 
##                                          2 
##                              Myristic acid 
##                                          2 
##                   gamma-Glutamylisoleucine 
##                                          9 
##              alpha-N-Phenylacetylglutamine 
##                                          1 
##                    Octadecatetraenoic acid 
##                                          2 
##                            Eicosenoic acid 
##                                          1 
##                      Docosapentaenoic acid 
##                                          8 
##                        Docosatrienoic acid 
##                                          2 
##                         Tetracosenoic acid 
##                                          2 
## 5alpha-Androstane-3alpha-ol-17-one sulfate 
##                                         18 
##                       Ursodeoxycholic acid 
##                                          6 
##                       Pregnenolone sulfate 
##                                          8 
##                                Cholic acid 
##                                         17 
## [1] "badi"
##    glutamine-13C5 [iSTD]          Lignoceric acid     Ursodeoxycholic acid 
##                        1                        2                        6 
## Taurohyodeoxycholic acid 
##                        1 
## [1] "wadi"
##         alanine-13C3 [iSTD]                 Ketoleucine 
##                           5                           3 
##        Hydroxyoctanoic acid              Kynurenic acid 
##                           1                           1 
##     citric acid-13C6 [iSTD]       N-Acetylphenylalanine 
##                           2                          14 
##      Hydroxydodecanoic acid   Hydroxytetradecanoic acid 
##                           2                           5 
##                  Xanthosine       Heptadecanedioic acid 
##                           1                           2 
##                Behenic acid Tauro-alpha-muricholic acid 
##                           1                           7 
##            Taurocholic acid 
##                           1

4 Imputation

4.1 Imputation (test samples)

pass1a.0.nr2 <- list()
for(i in 1:length(pass1a.0)){
  if(na_vals_impute[[i]] > 0){
    print(tis[[i]])
    glue("Features & Values to impute:")
    feature_impute[[i]]
    pass1a.0.nr2[[i]] <-impute.knn(pass1a.0.nr1[[i]],k=10)$data
    #view the features before and after imputation
    par(mfrow=c(2,1),bg="black")
    image(as.matrix(pass1a.0.nr1[[i]][,pass1a.0.nr1.f.c0[[i]]>0]),col=redblue100,axes=F)
    image(as.matrix(pass1a.0.nr2[[i]][, pass1a.0.nr1.f.c0[[i]]>0]),col=redblue100,axes=F)
    par(mfrow=c(1,1) ,bg="white")
    glue("Verified no missing values: {sum(is.na(pass1a.0.nr2[[i]]))}")  #verified 0
  }else{
    print('No missing values to impute')
    pass1a.0.nr2[[i]] <- pass1a.0.nr1[[i]]
  }
}

## [1] "No missing values to impute"
## [1] "hip"

## [1] "gas"

## [1] "hrt"

## [1] "kid"

## [1] "lun"

## [1] "liv"

## [1] "badi"

## [1] "wadi"

5 NxN heatmaps

5.1 Plotting NxN heatmaps (+ reference samples)

a <- list(0.82,0.85,0.80,
       0.80,0.80,0.82,
       0.80,0.80,0.80)
b <- 1
par(mfrow=c(3,3), mar = c(0,0,3,0))
i <- 1
for(i in 1:length(pass1a.0)){
  print(tis[[i]])
  x <- tmp.ref1a[[i]][,1]
  sidebar  <- ( (b - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  sidebar <- cbind(sidebar,sidebar,sidebar)
  cor.tmp<-cor(t(pass1a.0.1[[i]]),method="spearman",use="pairwise.complete.obs") #includes ref samples
  print(dim(cor.tmp))
  image(cbind(cor.tmp,sidebar),
            col=redblue100,axes=FALSE,zlim=c(a[[i]],1), main=glue("{TIS[[i]]}2-{NR1[[i]]},
                                                                  run-order, z=({a[[i]]},1)"), asp = 1)
}

## [1] "pla"
## [1] 96 96

## [1] "hip"
## [1] 113 113

## [1] "gas"
## [1] 98 98

## [1] "hrt"
## [1] 120 120

## [1] "kid"
## [1] 120 120

## [1] "lun"
## [1] 120 120

## [1] "liv"
## [1] 103 103

## [1] "badi"
## [1] 116 116

## [1] "wadi"
## [1] 112 112

# if(knit_time){
#   rpneg_pass1a.0.order[[i]] %>%
#     arrange(sample_order) %>%
#     select(sample_id, sample_type, sample_order) %>%
#     knitr::kable(format = "html") %>%
#     scroll_box(width = "100%", height = "400px")
# }

5.2 Plotting NxN heatmaps (test samples) and Identification of Outlier Samples

cor.tmp <- o.s <- list()
par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  x <- tmp.sample1a[[i]][,2]
  sidebar  <- ( (b - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  sidebar <- cbind(sidebar,sidebar,sidebar)
  cor.tmp[[i]]<-cor(t(pass1a.0.nr1[[i]]),method="spearman",use="pairwise.complete.obs")
  glue("Verified {N[[i]]} test samples: {dim(cor.tmp)[1]}") #verified, =78
  image(
    cbind(cor.tmp[[i]][order(tmp.sample1a[[i]][,2]),],sidebar),
    col=redblue100,axes=FALSE, zlim=c(a[[i]],1), main=glue("{TIS[[i]]}2-{N1[[i]]}, 
                                                           Run-Order, z=({a[[i]]},1)"), asp = 1)
}

# Outlier sample filter
o.f <- list(0.89, 0.915, 0.89,
            0.89, 0.90, 0.90,
            0.90, 0.89, 0.87)
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  plot(apply(cor.tmp[[i]][order(tmp.sample1a[[i]][,2]),order(tmp.sample1a[[i]][,2])],1,median), main = glue("{tis[[i]]}"),
       xlab = 'Samples (Run-Order)', ylab = 'Cor Medians'); abline(h=o.f[[i]], col = 'blue')
}

# Determine which samples are outliers
for(i in 1:length(pass1a.0)){
  if(is.na(o.f[[i]])){
    o.s[[i]] <- NA
  }else{
    o.s[[i]] <- colnames(cor.tmp[[i]])[apply(cor.tmp[[i]],1,median)<o.f[[i]]]
  }
  glue("Outlier Samples: {length(o.s[[i]])}")
}
# Examine outliers
glue("Outlier Samples in {TIS[[1]]}:")

## Outlier Samples in PLA:

rpneg_pass1a.0.order[[1]] %>% filter(sample_id %in% o.s[[1]]) %>%
  knitr::kable(format = "html") %>% scroll_box(width = "100%", height = "100%")

sample_id	sample_type	sample_order	Key.anirandgroup	Registration.sex
90010013104	Sample	31	Exercise - IPE	1

glue("Outlier Samples in {TIS[[2]]}:")

## Outlier Samples in HIP:

rpneg_pass1a.0.order[[2]] %>% filter(sample_id %in% o.s[[2]]) %>%
  knitr::kable(format = "html") %>% scroll_box(width = "100%", height = "100%")

sample_id	sample_type	sample_order	Key.anirandgroup	Registration.sex
90010015206	Sample	101	Exercise - IPE	1

glue("Outlier Samples in {TIS[[3]]}:")

## Outlier Samples in GAS:

rpneg_pass1a.0.order[[3]] %>% filter(sample_id %in% o.s[[3]]) %>%
  knitr::kable(format = "html") %>% scroll_box(width = "100%", height = "100%")

sample_id	sample_type	sample_order	Key.anirandgroup	Registration.sex
90136015506	Sample	25	Exercise - IPE	2
90156015506	Sample	53	Exercise - 1 hr	1
90007015506	Sample	57	Exercise - 0.5 hr	2
90010015506	Sample	63	Exercise - IPE	1
90012015506	Sample	77	Exercise - 1 hr	1
90025015506	Sample	85	Exercise - 0.5 hr	2

glue("Outlier Samples in {TIS[[4]]}:")

## Outlier Samples in HRT:

rpneg_pass1a.0.order[[4]] %>% filter(sample_id %in% o.s[[4]]) %>%
  knitr::kable(format = "html") %>% scroll_box(width = "100%", height = "100%")

sample_id	sample_type	sample_order	Key.anirandgroup	Registration.sex
90156015807	Sample	63	Exercise - 1 hr	1
90012015807	Sample	68	Exercise - 1 hr	1
90010015807	Sample	108	Exercise - IPE	1

glue("Outlier Samples in {TIS[[5]]}:")

## Outlier Samples in KID:

rpneg_pass1a.0.order[[5]] %>% filter(sample_id %in% o.s[[5]]) %>%
  knitr::kable(format = "html") %>% scroll_box(width = "100%", height = "100%")

sample_id	sample_type	sample_order	Key.anirandgroup	Registration.sex
90012015906	Sample	74	Exercise - 1 hr	1
90010015906	Sample	104	Exercise - IPE	1

glue("Outlier Samples in {TIS[[6]]}:")

## Outlier Samples in LUN:

rpneg_pass1a.0.order[[6]] %>% filter(sample_id %in% o.s[[6]]) %>%
  knitr::kable(format = "html") %>% scroll_box(width = "100%", height = "100%")

sample_id	sample_type	sample_order	Key.anirandgroup	Registration.sex
90010016606	Sample	69	Exercise - IPE	1
90012016606	Sample	98	Exercise - 1 hr	1

glue("Outlier Samples in {TIS[[7]]}:")

## Outlier Samples in LIV:

rpneg_pass1a.0.order[[7]] %>% filter(sample_id %in% o.s[[7]]) %>%
  knitr::kable(format = "html") %>% scroll_box(width = "100%", height = "100%")

sample_id	sample_type	sample_order	Key.anirandgroup	Registration.sex
90010016805	Sample	47	Exercise - IPE	1
90007016805	Sample	79	Exercise - 0.5 hr	2

glue("Outlier Samples in {TIS[[8]]}:")

## Outlier Samples in BADI:

rpneg_pass1a.0.order[[8]] %>% filter(sample_id %in% o.s[[8]]) %>%
  knitr::kable(format = "html") %>% scroll_box(width = "100%", height = "100%")

sample_id	sample_type	sample_order	Key.anirandgroup	Registration.sex
90012016906	Sample	31	Exercise - 1 hr	1
90010016906	Sample	33	Exercise - IPE	1
90156016906	Sample	72	Exercise - 1 hr	1

glue("Outlier Samples in {TIS[[9]]}:")

## Outlier Samples in WADI:

rpneg_pass1a.0.order[[9]] %>% filter(sample_id %in% o.s[[9]]) %>%
  knitr::kable(format = "html") %>% scroll_box(width = "100%", height = "100%")

sample_id	sample_type	sample_order	Key.anirandgroup	Registration.sex
90011017005	Sample	35	Exercise - 1 hr	2
90010017005	Sample	72	Exercise - IPE	1
90156017005	Sample	82	Exercise - 1 hr	1
90012017005	Sample	88	Exercise - 1 hr	1

5.3 Visualize the Sex-Time-Control NxN Heatmap

par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  x <- tmp.sample1a[[i]][,3]
  sex.type  <- ( (b - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  x <- tmp.sample1a[[i]][,4]
  control.type  <- ( (b - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  x <- as.factor(tmp.sample1a[[i]][,5]) %>% as.numeric()
  time.type  <- ( (b - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  sidebar<-cbind(sex.type,sex.type,sex.type, 
               control.type,control.type,control.type, 
               time.type,time.type,time.type)
  image(
    cbind(
    cor.tmp[[i]][order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5]),
          order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5])],
   sidebar[order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5]),]),
    col=redblue100,axes=FALSE, zlim=c(a[[i]],1), main=glue("{TIS[[i]]}2-{N[[i]]}, 
                                                           Sex-Control-Time, z=({a[[i]]},1)"), asp = 1)
}

par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  plot(apply(cor.tmp[[i]][order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5]), 
                         order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5])],1,median), main = glue("{tis[[i]]}"))
}

5.4 NxN Heatmap: Test and Reference Samples by Run Order

cor.tmp1 <- list()
par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  x <- is.na(tmp.ref1a[[i]][,2]) %>% as.integer()
  sample.type  <- ( (b - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  sidebar<-cbind(sample.type,sample.type,sample.type,sample.type,sample.type)
  cor.tmp1[[i]]<-cor(t(pass1a.0.1[[i]]),method="spearman",use="pairwise.complete.obs") #includes ref samples
  image(
    cbind(cor.tmp1[[i]], sidebar),
    col=redblue100,axes=FALSE, zlim=c(a[[i]],1), main=glue("{TIS[[i]]}2, 
                                                         Run-Order (+Ref-type), z=({a[[i]]},1)"), asp = 1)
}

5.5 Log Variables (3)

run_var <- list(0, 0, 0,
                0, 0, 0,
                0, 0, 0)
outlier_sample_n <- list( length(o.s[[1]]),length(o.s[[2]]),length(o.s[[3]]),
                         length(o.s[[4]]),length(o.s[[5]]),length(o.s[[6]]),
                         length(o.s[[7]]),length(o.s[[8]]),length(o.s[[9]]) )
outlier_samples <- o.s
outlier_filter <- o.f

6 N-P Heatmaps, Median and SD Distributions, & Transformations

6.0.1 N-P Heatmaps: Log2, imputed, sample order,

Note: Samples as columns Sidebar: p-values, t-values, and sig (binary) comparing normal test samples and outlier test samples

n.s <- hmo <- list()
par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  n.s[[i]] <- row.names(pass1a.0.nr2[[i]])[!row.names(pass1a.0.nr2[[i]]) %in% o.s[[i]]]
  hmo[[i]] <- heatmap(pass1a.0.nr2[[i]])$colInd
}

for(i in 1:length(pass1a.0)){
  image(pass1a.0.nr2[[i]][order(tmp.sample1a[[i]][,2]),hmo[[i]]]
      ,col=redblue100,axes=F,main=glue("{TIS[[i]]}2, 
                                       f{P1[[i]]}-HC, Run-Order"), asp = 1)
}

par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  plot(apply(pass1a.0.nr2[[i]][order(tmp.sample1a[[i]][,2]),hmo[[i]]],1,median), main = glue("{tis[[i]]}"))
}

6.0.2 N-P Heatmap: Log2, imputed, Sex-Control-Time order (2)

Note: pass1a.0.2 and pass1a.0.nr2 represent knn-imputed and log2 Note: Samples as columns

par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  a[[i]] <- range(pass1a.0.nr2[[i]])[1]
  b[[i]] <- range(pass1a.0.nr2[[i]])[2]
  x <- tmp.sample1a[[i]][,3]
  sex.type  <- ( (b[[i]] - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  x <- tmp.sample1a[[i]][,4]
  control.type <- ( (b[[i]] - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  x <- as.factor(tmp.sample1a[[i]][,5]) %>% as.numeric()
  time.type  <- ( (b[[i]] - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  sidebar<-cbind(sex.type,sex.type,sex.type,sex.type,sex.type, 
               control.type,control.type, control.type, control.type, control.type,
               time.type, time.type, time.type, time.type, time.type)
  image(
    cbind(pass1a.0.nr2[[i]][order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5]), hmo[[i]] ], 
      sidebar[order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5]),]),
      col=redblue100,axes=F,main=glue("{TIS[[i]]}2, 
                                    f{P1[[i]]}-HC, Sex-Control-Time"), asp =1)
}

par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  plot(apply(pass1a.0.nr2[[i]][order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5]),
     hmo[[i]]],1,median), main = glue("{tis[[i]]}"))
}

6.0.3 Transformations – only for the test samples. This can be revised to remove outlier samples

Strategy is not functionally programed (must revert log transformed back to linear for pass1a.0.3a)

pass1a.0.3a<-pass1a.0.3b<-pass1a.03c<-pass1a.03c2<-pass1a.02b<-pass1a.03d<-pass1a.03d2<-pass1a.02b2<-pass1a.03d3 <- list()
for(i in 1:length(pass1a.0)){
  # pass1a.0.r3a<-pass1a.0.r3b<-pass1a.0.r3c<-pass1a.0.r3c2<-pass1a.0.r2b<-pass1a.0.r3d<-pass1a.0.r3d2<-pass1a.0.r2b2<-pass1a.0.r3d3 <-pass1a.0.2
  pass1a.0.3a[[i]]<-pass1a.0.3b[[i]]<-pass1a.03c[[i]]<-pass1a.03c2[[i]]<-pass1a.02b[[i]]<-pass1a.03d[[i]]<-pass1a.03d2[[i]]<-pass1a.02b2[[i]]<-pass1a.03d3[[i]]<-pass1a.0.nr2[[i]]
  #pass1a.03c3<-pass1a.0.3b2<-pass1a.0.nr2
}

6.0.4 Examine the median and mean distributions (2)

par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  tmp.s.median <- apply(pass1a.0.nr2[[i]],1, median)
  tmp.s.mean <- apply(pass1a.0.nr2[[i]],1, mean)
  plot(tmp.s.median,tmp.s.mean, asp = 1, main = glue("{tis[[i]]}")); abline(0,1)
}

6.0.5 Examine the median and sd distribution (1)

tmp.f.median <- tmp.f.sd <- list()
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  tmp.f.median[[i]] <- apply(2^pass1a.0.nr2[[i]],2, median)
  tmp.f.sd[[i]] <- apply(2^pass1a.0.nr2[[i]],2, sd)
  plot(y = tmp.f.sd[[i]], x = tmp.f.median[[i]],log="xy", main = glue("{tis[[i]]}"))
}

6.0.6 Examine the median and sd distribution (1; w/wo outlier samples)

tmp2.f.median <- list()
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  tmp <- 2^pass1a.0.nr2[[i]][n.s[[i]],]
  #dim(tmp)
  tmp2.f.median[[i]] <- apply(tmp,2, median)
  #tmp2.f.sd <- apply(tmp,2, sd)
  plot(tmp.f.median[[i]],tmp2.f.median[[i]],log="xy", main = glue("{tis[[i]]}"))
  #plot(tmp.f.sd,tmp2.f.sd,log="xy")
}

tmp2.f.sd <- list()
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  tmp <- 2^pass1a.0.nr2[[i]][n.s[[i]],]
  tmp2.f.sd[[i]] <- apply(tmp,2, sd)
  plot(tmp.f.sd[[i]],tmp2.f.sd[[i]],log="xy",main = glue("{tis[[i]]}"))
}

6.1 Center by median, scale by standard deviation across features (3a)

Verification of the first feature

par(mfrow=c(3,3), mar = c(2,3,2,0))
j <- 148
for(i in 1:length(pass1a.0)){
  for (j in 1:length(tmp.f.median[[i]])) {
    pass1a.0.3a[[i]][,j]<-(2^pass1a.0.nr2[[i]][,j] - tmp.f.median[[i]][j])/ tmp.f.sd[[i]][j]
  }
  plot(2^pass1a.0.nr2[[i]][,1],pass1a.0.3a[[i]][,1], main = glue("{tis[[i]]}")) #spot-check, verified
}

6.1.1 N-P Heatmap: median-centered, sd-scaled imputed, run order (3a)

hmo2 <- list()
for(i in 1:length(pass1a.0)){
  hmo2[[i]] <- heatmap(pass1a.0.3a[[i]])$colInd
}

# Run Order; Original HC
###########################
par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  image(pass1a.0.3a[[i]][order(tmp.sample1a[[i]][,2]), hmo[[i]]],
    col=redblue100,axes=F,main=glue("{TIS[[i]]}3a, 
                                    f{P1[[i]]}-HC, Run-Order"), asp = 1)
}

# Run Order; New HC
###########################
par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  image(pass1a.0.3a[[i]][order(tmp.sample1a[[i]][,2]), hmo2[[i]]],
    col=redblue100,axes=F,main=glue("{TIS[[i]]}3a, 
                                    f{P1[[i]]}-HC(3a), Run-Order"), asp = 1)
}

par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  plot(apply(pass1a.0.3a[[i]][order(tmp.sample1a[[i]][,2]),
     hmo[[i]]],1,median), main = glue("{tis[[i]]}"), ylim = c(-2,2)); abline(h = 0, col = 'blue')
}

6.1.2 Examine the median and sd distribution (2)

tmp.f.median2 <- tmp.f.sd2 <- list()
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  tmp.f.median2[[i]] <- apply(pass1a.0.nr2[[i]],2, median)
  tmp.f.sd2[[i]] <- apply(pass1a.0.nr2[[i]],2, sd)
  plot(y = tmp.f.sd2[[i]], x = tmp.f.median2[[i]],log="xy", main = glue("{tis[[i]]}"))
}

6.1.3 Examine the median and sd distribution (2; w/wo outlier samples)

tmp2.f.median2 <- list()
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  tmp <-pass1a.0.nr2[[i]][n.s[[i]],]
  #dim(tmp)
  tmp2.f.median2[[i]] <- apply(tmp,2, median)
  plot(tmp.f.median2[[i]],tmp2.f.median2[[i]],log="xy", main = glue("{tis[[i]]}"))
}

tmp2.f.sd2 <- list()
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  tmp <- pass1a.0.nr2[[i]][n.s[[i]],]
  tmp2.f.sd2[[i]] <- apply(tmp,2, sd)
  plot(tmp.f.sd2[[i]],tmp2.f.sd2[[i]],log="xy",main = glue("{tis[[i]]}"))
}

6.2 Log2, Center by median, scale by standard deviation across features (3b)

Verification of the first & last feature

for(i in 1:length(pass1a.0)){
  for (j in 1:length(tmp2.f.median[[i]])) {
    pass1a.0.3b[[i]][,j]<-(pass1a.0.nr2[[i]][,j]- tmp.f.median2[[i]][j])/tmp.f.sd2[[i]][j]
  }
}
par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  plot(pass1a.0.nr2[[i]][,1],pass1a.0.3b[[i]][,1], main = glue("{tis[[i]]}")) #spot-check, verified
}

for(i in 1:length(pass1a.0)){
  plot(pass1a.0.nr2[[i]][,P1[[i]]],pass1a.0.3b[[i]][,P1[[i]]], main = glue("{tis[[i]]}")) #spot-check, verified
}

6.2.1 N-P Heatmap: median-centered, sd-scaled imputed, run order (3b)

hmo3 <- list()
for(i in 1:length(pass1a.0)){
  hmo3[[i]] <- heatmap(pass1a.0.3b[[i]])$colInd
}

# Run Order; Original HC
###########################
par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  image(pass1a.0.3b[[i]][order(tmp.sample1a[[i]][,2]), hmo[[i]]],
    col=redblue100,axes=F,main=glue("{TIS[[i]]}3b, 
                                    f{P1[[i]]}-HC, Run-Order"), asp = 1)
}

# Run Order; New HC
###########################
par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  image(pass1a.0.3b[[i]][order(tmp.sample1a[[i]][,2]), hmo3[[i]]],
    col=redblue100,axes=F,main=glue("{TIS[[i]]}3b, 
                                    f{P1[[i]]}-HC(3b), Run-Order"), asp = 1)
}

par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  plot(apply(pass1a.0.3b[[i]][order(tmp.sample1a[[i]][,2]),
     hmo[[i]]],1,median), main = glue("{tis[[i]]}"), ylim = c(-1.5,1.5)); abline(h = 0, col = 'blue')
}

6.2.2 N-P Heatmap: median-centered, sd-scaled imputed, Sex-Control-Time order (3b)

par(mfrow=c(3,3), mar = c(0,0,3,0))
for(i in 1:length(pass1a.0)){
  a[[i]] <- range(pass1a.0.3b[[i]])[1]
  b[[i]] <- range(pass1a.0.3b[[i]])[2]
  x <- tmp.sample1a[[i]][,3]
  sex.type  <- ( (b[[i]] - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  x <- tmp.sample1a[[i]][,4]
  control.type <- ( (b[[i]] - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  x <- as.factor(tmp.sample1a[[i]][,5]) %>% as.numeric()
  time.type  <- ( (b[[i]] - a[[i]]) * ((x - min(x))/((max(x) - min(x)))) + a[[i]]) 
  sidebar<-cbind(sex.type,sex.type,sex.type,sex.type,sex.type, 
               control.type,control.type, control.type, control.type, control.type,
               time.type, time.type, time.type, time.type, time.type)
  image(
    cbind(pass1a.0.3b[[i]][order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5]), hmo[[i]] ], 
      sidebar[order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5]),]),
      col=redblue100,axes=F,main=glue("{TIS[[i]]}2, 
                                    f{P1[[i]]}-HC, Sex-Control-Time"), asp =1)
}

par(mfrow=c(3,3), mar = c(2,3,2,0))
for(i in 1:length(pass1a.0)){
  plot(apply(pass1a.0.3b[[i]][order(tmp.sample1a[[i]][,3],tmp.sample1a[[i]][,4],tmp.sample1a[[i]][,5]),
     hmo[[i]]],1,median), main = glue("{tis[[i]]}"), ylim = c(-1.5,1.5)); abline(h = 0, col = 'blue')
}

6.2.3 Examine the sample median and sd distributions (3b)

par(mfrow=c(3,3), mar = c(0,2,1,0))
for(i in 1:length(pass1a.0)){
  boxplot(data.frame(t(pass1a.0.3b[[i]])), ylim = c(-2,2), main = glue("{tis[[i]]}"), xaxt = "n")
}

glue("Outlier samples removed:")

## Outlier samples removed:

par(mfrow=c(3,3), mar = c(0,2,1,0))
for(i in 1:length(pass1a.0)){
  boxplot(data.frame(t(pass1a.0.3b[[i]][n.s[[i]],])), ylim = c(-2,2), main = glue("{tis[[i]]}"), xaxt = "n")
}

6.2.4 TODO (in additional script): Collect sample medians independent of run-order outliers to determine if samples should be centered and/or scaled

6.2.5 TODO: Incorporate outlier sample removal

6.2.6 TODO: Incorporate flagged samples

6.2.7 TODO: Incorporate outlier feature removal from transformations

7 Save the processed data

7.1 Concatenate the Processing Decisions

pass1a.0.vars <- list()
comments <- list(NA, NA, NA,
                 NA, NA, NA,
                 NA, NA, NA)
pass1a.0.vars <- data.frame()
for(i in 1:length(pass1a.0)){
  feature_impute[[i]] = names(feature_impute[[i]])
  # The processing decisions
  pass1a.0.vars <- rbind(pass1a.0.vars, data.frame(ds, site, tech, 'tis' = tis[[i]], 'NR' = NR[[i]], 'N' = N[[i]], 'P' = P[[i]], 
                                   neg_vals, zero_vals, feature_na_filter, 'P1' = P1[[i]], 'NR1' = NR1[[i]], 
                                   'N1' = N1[[i]], 'imputed_features' = paste(feature_impute[[i]], collapse = '; '), 
                                   'na_vals_impute' = na_vals_impute[[i]], knn_k, 'run_var' = run_var[[i]], 
                                   'outlier_sample_n' = outlier_sample_n[[i]], 'outlier_samples' =  paste(o.s[[i]], collapse = '; '), 
                                   'outlier_filter' = o.f[[i]],  'internal_standards_n' = length(is[[i]]), 
                                   'internal_standards' = paste(is[[i]], collapse = '; '), 'comments' = comments[[i]]))
}

# visualize the processing decisions
if(knit_time){
  pass1a.0.vars %>%
    knitr::kable(format = "html") %>%
    scroll_box(width = "100%", height = "100%")
}

ds	site	tech	tis	NR	N	P	feature_na_filter	P1	NR1	N1	imputed_features	na_vals_impute	knn_k	outlier_sample_n	outlier_samples	outlier_filter	internal_standards_n	internal_standards	comments
pass1a	umich	rpneg	pla	96	77	97	20	96	96	77		0	10	1	90010013104	0.890	10	valine-13C5 [iSTD]; aspartate-13C4 [iSTD]; glutamine-13C5 [iSTD]; glutamate-13C5 [iSTD]; phenylalanine-13C9 [iSTD]; tyrosine-13C9 [iSTD]; citric acid-13C6 [iSTD]; tryptophan-15N2 [iSTD]; Gibberellic acid [iSTD]; 24-epi-Brassinolide [iSTD]	NA
pass1a	umich	rpneg	hip	113	78	139	20	109	113	78	glutamate-13C5 [iSTD]	3	10	1	90010015206	0.915	22	alanine-13C3 [iSTD]; lactate-13C3 [iSTD]; serine-13C3 [iSTD]; valine-13C5 [iSTD]; threonine-13C4 [iSTD]; Asparagine-15N2 [iSTD]; aspartate-13C4 [iSTD]; isoleucine-13C6 [iSTD]; Leucine-13C/Isoleucine-13C6 [iSTD]; leucine-13C6 [iSTD]; malate-13C4 [iSTD]; anthranilic acid-15N [iSTD]; glutamine-13C5 [iSTD]; glutamate-13C5 [iSTD]; lysine-13C6 [iSTD]; histidine-13C6 [iSTD]; phenylalanine-13C9 [iSTD]; tyrosine-13C9 [iSTD]; citric acid-13C6 [iSTD]; tryptophan-15N2 [iSTD]; Gibberellic acid [iSTD]; 24-epi-Brassinolide [iSTD]	NA
pass1a	umich	rpneg	gas	98	78	78	20	76	98	78	Tetradecadienoic acid; Myristoleic acid; Deoxycytidine; Pentadecylic acid; Xanthosine	20	10	6	90136015506; 90156015506; 90007015506; 90010015506; 90012015506; 90025015506	0.890	0		NA
pass1a	umich	rpneg	hrt	120	78	94	20	94	120	78	alanine-13C3 [iSTD]; 3-Methyl-2-oxovaleric acid; Ketoleucine; Lauric acid; Lignoceric acid; Ursodeoxycholic acid	19	10	3	90156015807; 90012015807; 90010015807	0.890	12	alanine-13C3 [iSTD]; serine-13C3 [iSTD]; valine-13C5 [iSTD]; aspartate-13C4 [iSTD]; glutamine-13C5 [iSTD]; glutamate-13C5 [iSTD]; phenylalanine-13C9 [iSTD]; tyrosine-13C9 [iSTD]; citric acid-13C6 [iSTD]; tryptophan-15N2 [iSTD]; Gibberellic acid [iSTD]; 24-epi-Brassinolide [iSTD]	NA
pass1a	umich	rpneg	kid	120	78	87	20	87	120	78	Lactic acid; Capric acid; Pregnenolone sulfate; Cholic acid	9	10	2	90012015906; 90010015906	0.900	10	alanine-13C3 [iSTD]; valine-13C5 [iSTD]; aspartate-13C4 [iSTD]; glutamine-13C5 [iSTD]; glutamate-13C5 [iSTD]; tyrosine-13C9 [iSTD]; citric acid-13C6 [iSTD]; tryptophan-15N2 [iSTD]; Gibberellic acid [iSTD]; 24-epi-Brassinolide [iSTD]	NA
pass1a	umich	rpneg	lun	120	78	102	20	102	120	78	Tauro-alpha-muricholic acid	5	10	2	90010016606; 90012016606	0.900	16	alanine-13C3 [iSTD]; lactate-13C3 [iSTD]; aspartate-13C4 [iSTD]; isoleucine-13C6 [iSTD]; leucine-13C6 [iSTD]; anthranilic acid-15N [iSTD]; glutamine-13C5 [iSTD]; glutamate-13C5 [iSTD]; lysine-13C6 [iSTD]; histidine-13C6 [iSTD]; phenylalanine-13C9 [iSTD]; tyrosine-13C9 [iSTD]; citric acid-13C6 [iSTD]; tryptophan-15N2 [iSTD]; Gibberellic acid [iSTD]; 24-epi-Brassinolide [iSTD]	NA
pass1a	umich	rpneg	liv	103	78	81	20	81	103	78	Tetradecadienoic acid; Myristoleic acid; Myristic acid; gamma-Glutamylisoleucine; alpha-N-Phenylacetylglutamine; Octadecatetraenoic acid; Eicosenoic acid; Docosapentaenoic acid; Docosatrienoic acid; Tetracosenoic acid; 5alpha-Androstane-3alpha-ol-17-one sulfate; Ursodeoxycholic acid; Pregnenolone sulfate; Cholic acid	80	10	2	90010016805; 90007016805	0.900	0		NA
pass1a	umich	rpneg	badi	116	78	143	20	123	116	78	glutamine-13C5 [iSTD]; Lignoceric acid; Ursodeoxycholic acid; Taurohyodeoxycholic acid	10	10	3	90012016906; 90010016906; 90156016906	0.890	22	alanine-13C3 [iSTD]; lactate-13C3 [iSTD]; serine-13C3 [iSTD]; valine-13C5 [iSTD]; threonine-13C4 [iSTD]; Asparagine-15N2 [iSTD]; aspartate-13C4 [iSTD]; isoleucine-13C6 [iSTD]; Leucine-13C/Isoleucine-13C6 [iSTD]; leucine-13C6 [iSTD]; malate-13C4 [iSTD]; anthranilic acid-15N [iSTD]; glutamine-13C5 [iSTD]; glutamate-13C5 [iSTD]; lysine-13C6 [iSTD]; histidine-13C6 [iSTD]; phenylalanine-13C9 [iSTD]; tyrosine-13C9 [iSTD]; citric acid-13C6 [iSTD]; tryptophan-15N2 [iSTD]; Gibberellic acid [iSTD]; 24-epi-Brassinolide [iSTD]	NA
pass1a	umich	rpneg	wadi	112	76	78	20	76	112	76	alanine-13C3 [iSTD]; Ketoleucine; Hydroxyoctanoic acid; Kynurenic acid; citric acid-13C6 [iSTD]; N-Acetylphenylalanine; Hydroxydodecanoic acid; Hydroxytetradecanoic acid; Xanthosine; Heptadecanedioic acid; Behenic acid; Tauro-alpha-muricholic acid; Taurocholic acid	45	10	4	90011017005; 90010017005; 90156017005; 90012017005	0.870	9	alanine-13C3 [iSTD]; valine-13C5 [iSTD]; aspartate-13C4 [iSTD]; glutamine-13C5 [iSTD]; glutamate-13C5 [iSTD]; citric acid-13C6 [iSTD]; tryptophan-15N2 [iSTD]; Gibberellic acid [iSTD]; 24-epi-Brassinolide [iSTD]	NA

7.2 Save the Data

# pla
pla1a.0.nr1 <- pass1a.0.nr1[[1]]
pla1a.0.1 <- pass1a.0.1[[1]]
pla1a.0.nr2 <- pass1a.0.nr2[[1]]
pla1a.0.3a <- pass1a.0.3a[[1]]
pla1a.0.3b <- pass1a.0.3b[[1]]
pla1a.0.vars <- pass1a.0.vars[1,]
# save the data
save(pla1a.0.nr1, pla1a.0.1, pla1a.0.nr2, pla1a.0.3a, pla1a.0.3b, pla1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_pla1a.0.RData"))

# hip
hip1a.0.nr1 <- pass1a.0.nr1[[2]]
hip1a.0.1 <- pass1a.0.1[[2]]
hip1a.0.nr2 <- pass1a.0.nr2[[2]]
hip1a.0.3a <- pass1a.0.3a[[2]]
hip1a.0.3b <- pass1a.0.3b[[2]]
hip1a.0.vars <- pass1a.0.vars[2,]
# save the data
save(hip1a.0.nr1, hip1a.0.1, hip1a.0.nr2, hip1a.0.3a, hip1a.0.3b, hip1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_hip1a.0.RData"))

# gas
gas1a.0.nr1 <- pass1a.0.nr1[[3]]
gas1a.0.1 <- pass1a.0.1[[3]]
gas1a.0.nr2 <- pass1a.0.nr2[[3]]
gas1a.0.3a <- pass1a.0.3a[[3]]
gas1a.0.3b <- pass1a.0.3b[[3]]
gas1a.0.vars <- pass1a.0.vars[3,]
# save the data
save(gas1a.0.nr1, gas1a.0.1, gas1a.0.nr2, gas1a.0.3a, gas1a.0.3b, gas1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_gas1a.0.RData"))

# hrt
hrt1a.0.nr1 <- pass1a.0.nr1[[4]]
hrt1a.0.1 <- pass1a.0.1[[4]]
hrt1a.0.nr2 <- pass1a.0.nr2[[4]]
hrt1a.0.3a <- pass1a.0.3a[[4]]
hrt1a.0.3b <- pass1a.0.3b[[4]]
hrt1a.0.vars <- pass1a.0.vars[4,]
# save the data
save(hrt1a.0.nr1, hrt1a.0.1, hrt1a.0.nr2, hrt1a.0.3a, hrt1a.0.3b, hrt1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_hrt1a.0.RData"))

# kid
kid1a.0.nr1 <- pass1a.0.nr1[[5]]
kid1a.0.1 <- pass1a.0.1[[5]]
kid1a.0.nr2 <- pass1a.0.nr2[[5]]
kid1a.0.3a <- pass1a.0.3a[[5]]
kid1a.0.3b <- pass1a.0.3b[[5]]
kid1a.0.vars <- pass1a.0.vars[5,]
# save the data
save(kid1a.0.nr1, kid1a.0.1, kid1a.0.nr2, kid1a.0.3a, kid1a.0.3b, kid1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_kid1a.0.RData"))

# lun
lun1a.0.nr1 <- pass1a.0.nr1[[6]]
lun1a.0.1 <- pass1a.0.1[[6]]
lun1a.0.nr2 <- pass1a.0.nr2[[6]]
lun1a.0.3a <- pass1a.0.3a[[6]]
lun1a.0.3b <- pass1a.0.3b[[6]]
lun1a.0.vars <- pass1a.0.vars[6,]
# save the data
save(lun1a.0.nr1, lun1a.0.1, lun1a.0.nr2, lun1a.0.3a, lun1a.0.3b, lun1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_lun1a.0.RData"))

# liv
liv1a.0.nr1 <- pass1a.0.nr1[[7]]
liv1a.0.1 <- pass1a.0.1[[7]]
liv1a.0.nr2 <- pass1a.0.nr2[[7]]
liv1a.0.3a <- pass1a.0.3a[[7]]
liv1a.0.3b <- pass1a.0.3b[[7]]
liv1a.0.vars <- pass1a.0.vars[7,]
# save the data
save(liv1a.0.nr1, liv1a.0.1, liv1a.0.nr2, liv1a.0.3a, liv1a.0.3b, liv1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_liv1a.0.RData"))

# badi
badi1a.0.nr1 <- pass1a.0.nr1[[8]]
badi1a.0.1 <- pass1a.0.1[[8]]
badi1a.0.nr2 <- pass1a.0.nr2[[8]]
badi1a.0.3a <- pass1a.0.3a[[8]]
badi1a.0.3b <- pass1a.0.3b[[8]]
badi1a.0.vars <- pass1a.0.vars[8,]
# save the data
save(badi1a.0.nr1, badi1a.0.1, badi1a.0.nr2, badi1a.0.3a, badi1a.0.3b, badi1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_badi1a.0.RData"))

# wadi
wadi1a.0.nr1 <- pass1a.0.nr1[[9]]
wadi1a.0.1 <- pass1a.0.1[[9]]
wadi1a.0.nr2 <- pass1a.0.nr2[[9]]
wadi1a.0.3a <- pass1a.0.3a[[9]]
wadi1a.0.3b <- pass1a.0.3b[[9]]
wadi1a.0.vars <- pass1a.0.vars[9,]
# save the data
save(wadi1a.0.nr1, wadi1a.0.1, wadi1a.0.nr2, wadi1a.0.3a, wadi1a.0.3b, wadi1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_wadi1a.0.RData"))

# all tissues
save(pla1a.0.nr1, pla1a.0.1, pla1a.0.nr2, pla1a.0.3a, pla1a.0.3b, pla1a.0.vars,
     hip1a.0.nr1, hip1a.0.1, hip1a.0.nr2, hip1a.0.3a, hip1a.0.3b, hip1a.0.vars,
     gas1a.0.nr1, gas1a.0.1, gas1a.0.nr2, gas1a.0.3a, gas1a.0.3b, gas1a.0.vars,
     hrt1a.0.nr1, hrt1a.0.1, hrt1a.0.nr2, hrt1a.0.3a, hrt1a.0.3b, hrt1a.0.vars,
     kid1a.0.nr1, kid1a.0.1, kid1a.0.nr2, kid1a.0.3a, kid1a.0.3b, kid1a.0.vars,
     lun1a.0.nr1, lun1a.0.1, lun1a.0.nr2, lun1a.0.3a, lun1a.0.3b, lun1a.0.vars,
     liv1a.0.nr1, liv1a.0.1, liv1a.0.nr2, liv1a.0.3a, liv1a.0.3b, liv1a.0.vars,
     badi1a.0.nr1, badi1a.0.1, badi1a.0.nr2, badi1a.0.3a, badi1a.0.3b, badi1a.0.vars,
     wadi1a.0.nr1, wadi1a.0.1, wadi1a.0.nr2, wadi1a.0.3a, wadi1a.0.3b, wadi1a.0.vars,
  file = glue("{WD}/data/UM-rpneg/UM_rpneg_processed_pass1a.0.RData"))

8 Session Info

warnings()
session_info()

## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 4.0.2 (2020-06-22)
##  os       macOS  10.16                
##  system   x86_64, darwin17.0          
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       America/Toronto             
##  date     2021-12-15                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package      * version  date       lib source                                
##  assertthat     0.2.1    2019-03-21 [2] CRAN (R 4.0.2)                        
##  backports      1.2.1    2020-12-09 [2] CRAN (R 4.0.2)                        
##  broom          0.7.8    2021-06-24 [1] CRAN (R 4.0.2)                        
##  cachem         1.0.3    2021-02-04 [2] CRAN (R 4.0.2)                        
##  callr          3.7.0    2021-04-20 [1] CRAN (R 4.0.2)                        
##  cellranger     1.1.0    2016-07-27 [2] CRAN (R 4.0.2)                        
##  cli            3.0.1    2021-07-17 [1] CRAN (R 4.0.2)                        
##  coda           0.19-4   2020-09-30 [1] CRAN (R 4.0.2)                        
##  codetools      0.2-18   2020-11-04 [2] CRAN (R 4.0.2)                        
##  colorspace     2.0-2    2021-06-24 [1] CRAN (R 4.0.2)                        
##  crayon         1.4.1    2021-02-08 [2] CRAN (R 4.0.2)                        
##  curl           4.3.2    2021-06-23 [1] CRAN (R 4.0.2)                        
##  data.table     1.14.0   2021-02-21 [1] CRAN (R 4.0.2)                        
##  DBI            1.1.1    2021-01-15 [2] CRAN (R 4.0.2)                        
##  dbplyr         2.1.1    2021-04-06 [1] CRAN (R 4.0.2)                        
##  desc           1.3.0    2021-03-05 [1] CRAN (R 4.0.2)                        
##  devtools     * 2.4.2    2021-06-07 [1] CRAN (R 4.0.2)                        
##  digest         0.6.27   2020-10-24 [2] CRAN (R 4.0.2)                        
##  dplyr        * 1.0.7    2021-06-18 [1] CRAN (R 4.0.2)                        
##  ellipsis       0.3.2    2021-04-29 [1] CRAN (R 4.0.2)                        
##  evaluate       0.14     2019-05-28 [2] CRAN (R 4.0.1)                        
##  fansi          0.5.0    2021-05-25 [1] CRAN (R 4.0.2)                        
##  fastmap        1.1.0    2021-01-25 [2] CRAN (R 4.0.2)                        
##  forcats      * 0.5.1    2021-01-27 [2] CRAN (R 4.0.2)                        
##  fs             1.5.0    2020-07-31 [2] CRAN (R 4.0.2)                        
##  generics       0.1.0    2020-10-31 [2] CRAN (R 4.0.2)                        
##  ggfittext      0.9.1    2021-01-30 [2] CRAN (R 4.0.2)                        
##  ggplot2      * 3.3.5    2021-06-25 [1] CRAN (R 4.0.2)                        
##  glue         * 1.4.2    2020-08-27 [2] CRAN (R 4.0.2)                        
##  gridExtra      2.3      2017-09-09 [2] CRAN (R 4.0.2)                        
##  gtable         0.3.0    2019-03-25 [2] CRAN (R 4.0.2)                        
##  haven          2.3.1    2020-06-01 [2] CRAN (R 4.0.2)                        
##  highr          0.9      2021-04-16 [1] CRAN (R 4.0.2)                        
##  hms            1.1.0    2021-05-17 [1] CRAN (R 4.0.2)                        
##  htmltools      0.5.1.1  2021-01-22 [2] CRAN (R 4.0.2)                        
##  httr           1.4.2    2020-07-20 [2] CRAN (R 4.0.2)                        
##  impute       * 1.62.0   2020-04-27 [1] Bioconductor                          
##  inline         0.3.19   2021-05-31 [1] CRAN (R 4.0.2)                        
##  inspectdf      0.0.11   2021-04-02 [1] CRAN (R 4.0.2)                        
##  jsonlite       1.7.2    2020-12-09 [2] CRAN (R 4.0.2)                        
##  kableExtra   * 1.3.1    2020-10-22 [2] CRAN (R 4.0.2)                        
##  knitr          1.33     2021-04-24 [1] CRAN (R 4.0.2)                        
##  lattice        0.20-41  2020-04-02 [2] CRAN (R 4.0.2)                        
##  lifecycle      1.0.0    2021-02-15 [1] CRAN (R 4.0.2)                        
##  loo            2.4.1    2020-12-09 [1] CRAN (R 4.0.2)                        
##  lubridate      1.7.10   2021-02-26 [1] CRAN (R 4.0.2)                        
##  magrittr       2.0.1    2020-11-17 [2] CRAN (R 4.0.2)                        
##  MASS           7.3-53   2020-09-09 [2] CRAN (R 4.0.2)                        
##  matrixStats    0.60.0   2021-07-26 [1] CRAN (R 4.0.2)                        
##  memoise        2.0.0    2021-01-26 [2] CRAN (R 4.0.2)                        
##  modelr         0.1.8    2020-05-19 [2] CRAN (R 4.0.2)                        
##  MotrpacBicQC * 0.5.2    2021-06-25 [1] Github (MoTrPAC/MotrpacBicQC@43fb293) 
##  munsell        0.5.0    2018-06-12 [2] CRAN (R 4.0.2)                        
##  mvtnorm        1.1-2    2021-06-07 [1] CRAN (R 4.0.2)                        
##  naniar         0.6.1    2021-05-14 [1] CRAN (R 4.0.2)                        
##  pillar         1.6.2    2021-07-29 [1] CRAN (R 4.0.2)                        
##  pkgbuild       1.2.0    2020-12-15 [2] CRAN (R 4.0.2)                        
##  pkgconfig      2.0.3    2019-09-22 [2] CRAN (R 4.0.2)                        
##  pkgload        1.2.1    2021-04-06 [1] CRAN (R 4.0.2)                        
##  plyr           1.8.6    2020-03-03 [2] CRAN (R 4.0.2)                        
##  prettyunits    1.1.1    2020-01-24 [2] CRAN (R 4.0.2)                        
##  processx       3.5.2    2021-04-30 [1] CRAN (R 4.0.2)                        
##  progress       1.2.2    2019-05-16 [2] CRAN (R 4.0.2)                        
##  ps             1.6.0    2021-02-28 [1] CRAN (R 4.0.2)                        
##  purrr        * 0.3.4    2020-04-17 [2] CRAN (R 4.0.2)                        
##  R6             2.5.0    2020-10-28 [2] CRAN (R 4.0.2)                        
##  Rcpp           1.0.7    2021-07-07 [1] CRAN (R 4.0.2)                        
##  RcppParallel   5.1.4    2021-05-04 [1] CRAN (R 4.0.2)                        
##  readr        * 1.4.0    2020-10-05 [2] CRAN (R 4.0.2)                        
##  readxl         1.3.1    2019-03-13 [2] CRAN (R 4.0.2)                        
##  remotes        2.4.0    2021-06-02 [1] CRAN (R 4.0.2)                        
##  reprex         2.0.0    2021-04-02 [1] CRAN (R 4.0.2)                        
##  reshape2       1.4.4    2020-04-09 [2] CRAN (R 4.0.2)                        
##  rethinking   * 2.13     2021-08-13 [1] Github (rmcelreath/rethinking@2acf2fd)
##  rlang          0.4.11   2021-04-30 [1] CRAN (R 4.0.2)                        
##  rmarkdown      2.6      2020-12-14 [2] CRAN (R 4.0.2)                        
##  rprojroot      2.0.2    2020-11-15 [2] CRAN (R 4.0.2)                        
##  rstan        * 2.21.2   2020-07-27 [1] CRAN (R 4.0.2)                        
##  rstudioapi     0.13     2020-11-12 [2] CRAN (R 4.0.2)                        
##  rvest          1.0.0    2021-03-09 [1] CRAN (R 4.0.2)                        
##  scales         1.1.1    2020-05-11 [2] CRAN (R 4.0.2)                        
##  sessioninfo    1.1.1    2018-11-05 [2] CRAN (R 4.0.2)                        
##  shape          1.4.6    2021-05-19 [1] CRAN (R 4.0.2)                        
##  StanHeaders  * 2.21.0-7 2020-12-17 [1] CRAN (R 4.0.2)                        
##  stringi        1.7.3    2021-07-16 [1] CRAN (R 4.0.2)                        
##  stringr      * 1.4.0    2019-02-10 [2] CRAN (R 4.0.2)                        
##  testthat       3.0.4    2021-07-01 [1] CRAN (R 4.0.2)                        
##  tibble       * 3.1.3    2021-07-23 [1] CRAN (R 4.0.2)                        
##  tidyr        * 1.1.3    2021-03-03 [1] CRAN (R 4.0.2)                        
##  tidyselect     1.1.1    2021-04-30 [1] CRAN (R 4.0.2)                        
##  tidyverse    * 1.3.1    2021-04-15 [1] CRAN (R 4.0.2)                        
##  usethis      * 2.0.1    2021-02-10 [1] CRAN (R 4.0.2)                        
##  utf8           1.2.2    2021-07-24 [1] CRAN (R 4.0.2)                        
##  V8             3.4.2    2021-05-01 [1] CRAN (R 4.0.2)                        
##  vctrs          0.3.8    2021-04-29 [1] CRAN (R 4.0.2)                        
##  viridisLite    0.4.0    2021-04-13 [1] CRAN (R 4.0.2)                        
##  visdat         0.5.3    2019-02-15 [2] CRAN (R 4.0.2)                        
##  webshot        0.5.2    2019-11-22 [2] CRAN (R 4.0.2)                        
##  withr          2.4.2    2021-04-18 [1] CRAN (R 4.0.2)                        
##  xfun           0.24     2021-06-15 [1] CRAN (R 4.0.2)                        
##  xml2           1.3.2    2020-04-23 [2] CRAN (R 4.0.2)                        
##  yaml           2.2.1    2020-02-01 [2] CRAN (R 4.0.2)                        
## 
## [1] /Users/Alec/Library/R/4.0/library
## [2] /Library/Frameworks/R.framework/Versions/4.0/Resources/library

PASS1A/PASS1C Metabolomics: Data EDA/Cleaning/Imputation/Transformation–all_rpneg_pass1a.0 (pass1a 1.0 post-IS)

Alec Steep, Jun Li

12/12/2021